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We present a series of theoretical studies of the boundary between a superfluid and normal region 
in a partially polarized gas of strongly interacting fermions. We present mean-field estimates of the 
surface energy in this boundary as a function of temperature and scattering length. We discuss the 
structure of the domain wall, and use a previously introduced phenomonological model to study its 
influence on experimental observables. Our microscopic mean-field calculations are not consistent 
with the magnitude of the surface tension found from our phenomonological modelling of data from 
the Rice experiments. We conclude that one must search for novel mechanisms to explain the 
experiments. 



I. INTRODUCTION 

What happens when one tries to polarize a fermionic 
superfluid? Experiments at MIT and Rice have shown 
that when the fermions are interacting via resonant short 
range interactions, the fluid responds by phase separat- 
ing into a largely unpolarized superfluid region and a less 
polarized normal region p], 0, y, 0, [f| . The Rice exper- 
iments [3, HJ show a dramatic distortion of the central 
superfluid region in their trapped gas, pointing to signif- 
icant surface tension in the boundary. Here we present a 
theoretical study of this boundary and discuss the con- 
sequences of surface tension. 

The phase separation seen in these experiments arises 
because a zero temperature conventional s-wavc super- 
fluid is unable to accommodate spin polarization: all of 
the atoms in one spin state {]) are paired with atoms 
of the opposite spin {[). Changing the density ratio 
n^/ni from unity requires adding sufficient energy to 
break these pairs. Consequently, when excess particles 
of one spin state are added to a paired atomic cloud, 
those particles simply float to the surface, forming a nor- 
mal fluid. Given that there is a sharp boundary between 
the superfluid and normal region, the order parameter 
must vary rapidly, producing a surface energy. In sec- 
tion|TT]we use the Bogoliubov-de-Gennes(BdG) equations 
as a microscopic theory of this interface, and extract a 
dimcnsionlcss measure of the surface tension rj. 

In addition to the phase separation scenario seen in ex- 
periments, which was essentially predicted by Clogston 
and Chandrasekhar 0, 0], there have been many theo- 
retical proposals for how the superfluid can accomodate 
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HO, Hil [H, miHHi, Hi- We wiU make some com- 
ments on one class of states, the "FFLO" states intro- 
duced by Fulde, Ferrell, Larkin and Ovchinnikov d [t| 
where polarization resides in a periodic array of nodes 
in the superfluid order parameter, or is accommodated 
by creating supercurrents. Our mean field approach is 
sufficiently general that if such a phase existed in the 
parameter regime that we discuss, we could observe it. 



Like the experiments, and previous mean-field calcula- 
tions [281 ] . we see no sign of this phase near unitarity 
(where the scattering length is infinite). Interestingly, a 
recent density functional calculation [23 ] raises the possi- 
bility that the mean-field theory may be underestimating 
the stability of the FFLO state. We anticipate that in 
the near future more sophisticated numerical techniques 
will be able to unambiguously address the presence of the 
FFLO state. 

Like other theoretical calculations base d up on the Bo- 
goliubov de Gcnnes equations [13, HH, HI HJIH, I^E I^E 
l37l |38| we find that the order parameter and polarization 
oscillate in the domain wall separating the two regions. 
These oscillations, which in no way should be interpreted 
as an intervening phase [36| . are small at unitarity, but 
become larger as one approaches the BCS limit (small 
and negative scattering length a). At sufficiently small 
negative scattering length the decay length of these os- 
cillations diverge, signaling the onset of the bulk FFLO 
phase. The topology of the phase diagram of polar- 
ized fermions [22j, |39( at zero temperature features tri- 
critical points in both the BEC (1/kfa S> 1) and BCS 
(—1/kfCL 3> 1) limit, where the first order phase transi- 
tion between superfluid and normal state turn second or- 
der. As surface tension vanishes at the tricritical points, 
it reaches a maximum in the crossover region. 

In section IIIII we explore the consequences of surface 
tension on the shape of the superfluid-normal boundary 
for a unitary gas in an anisotropic harmonic trap. As in 
previous work [H El, EJ we use a simple elastic model 
for the boundary. By expanding the shape of the bound- 
ary in a Fourier series, we are able to compare its detailed 
structure with that of experiments. While our model ap- 
pears to capture a great many of the experimental fea- 
tures, it yields sharper density features. We attribute 
the discrepancies to the fact that we model the trapping 
potential as harmonic. 

This method complements the approach of Natu and 
Mueller [42j where the conditions of hydrostatic equilib- 
rium were used to derive a differential equation for the 
boundary. Like Haque et al. [4l[, we find that as one 
increases the surface tension, the boundary distorts from 
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an ellipse into a "capsule-like" shape. The most relevant 
experiments, performed at Rice in 2006 0, use small 
numbers of particles in a highly anisotropic trap. Due to 
the large surface area to volume ratio these experiments 
observe large distortions, consistent with surface tension 
which is one order of magnitude larger than predicted by 
our microscopic arguments. This behavior should be con- 
trasted with experimental studies at MIT Q which find 
no observable distortion of the superfluid-normal bound- 
ary. Taking into account the much smaller surface area 
to volume ratio in these experiments, this null observa- 
tion bounds 77 to be not much larger than the value we 
calculate. We have no explanation for this apparent dis- 
crepency. It is undoubtedly related to the fact that the 
Rice experiment finds a normal fluid whose local polariza- 
tion is almost 100%, while the normal fluid seen at MIT 
is always partially polarized, even at the lowest temper- 
atures. 



A. Background 

In the experiments of interest, fermionic alkali atoms 
(typically 6 Li) are trapped in a nominally harmonic op- 
tical potential. The atoms are transfered into two colli- 
sionally stable hyperfine states so that no spin relaxation 
occurs on the timescale of the experiment: the number 
of particles in each of the two spin states are separately 
conserved. Due to the short range nature of the interac- 
tions, at sufficiently low temperature scattering is forbid- 
den between two fermions in the same spin state. Hence 
interactions can be parameterized by a single scattering 
leng th a, which is a function of the applied magnetic field 

Sail. 

For a < the low energy scattering is attractive, and 
at low temperature the spin balanced system is a BCS 
supcrfluid. For a > the low energy scattering is re- 
pulsive, however there exists a two-body bound state in 
vacuum. The low temperature phase in this case is a 
Bose condensate of pairs. One of the remarkable predic- 
tions born out by experiments is that these two super- 
fluid pha ses are continuously connected to one another 
|, 111, |H, EE El, HI HI- Most interest has focused 
around the unitary point (a = ±00) where the scattering 
cross-section is maximal and in free space there exists 
a bound state at threshold. At this point the interac- 
tions do not provide a length-scale to the system, leading 
to universal thermodynamics where all thermodynamic 
functions can be expressed in terms of a power of the 
density times a universal function of the density scaled 
by the thermal wavel eng th nA 3 and the spin imbalance 
n t Am (EH, lED, HH, HH l55l | . In particular, if in the absence 
of a trap there is a flat phase boundary between a super- 
fluid and a normal region, with equ al p ressures on each 
side of the boundary, we showed in [4fJ] that any surface 
tension can be written as 

- = vf</3 (1) 
2m 



where n s is the density on the superfluid side of the 
boundary. A completely equivalent parameterization was 
later used by Haque and Stoof , 

To convert between the two parameterizations one needs 
to know the equations of state of the superfluid - which 
on dimensional grounds has the simple form 

where fi is the average chemical potential and is a 
dimensionless universal many body parameter [Fll . IH2I 
[HI, [55[ . According to quantum Monte-Carlo calculations 
P « -0.58 [H, HI, Ml, Hl| , which gives r\ = 8.1t7 s . 

If, as in the experiments, the boundary is curved, one 
expects to observe a pressure drop related to the curva- 
ture [42| . In that case, the dimensional argument leading 
to Eq. ([1]) yields an extra parameter, and one must take 
r/ = rj(Ap/p) to be a function of the relative pressure 
drop. As in all previous treatments, we will neglect the 
Ap dependance of this parameter. 

The presence of a superfluid-normal phase boundary 
for the trapped gas is understood by examining the fi- 
nite temperature phase diagram of a uniform Fermi gas, 
as sketched in figure[IJa) in the polarization-temperature 
plane. At zero temperature the superfluid and normal 
state are separated by a polarization driven first order 
phase transition. Since the polarization changes discon- 
tinuously, there is a "forbidden region" analogous to the 
one occuring in the density-temperature phase diagram 
for a liquid-gas phase transition. The location of the 
phase boundaries are found by a standard Maxwell con- 
struction, where the pressures (free energies) of the two 
phases are equated. When placed in a fixed volume (or 
confined by a harmonic trap) this first order phase tran- 
sition leads to a regime of phase coexistence, just as in 
the more familiar situation of a liquid-vapor transition. 

This first order phase transition is only found at suf- 
ficiently low temperatures. There is a tricritical point, 
above which the boundary becomes a continuous second 
order line. This behavior is consistent with the standard 
model of an unpolarized superfluid, which has a tempera- 
ture driven second order phase transition. This structure 
was experimentally mapped out in Q . 

In figure [ljc) the same diagram is shown at fixed 
temperature as a function of the chemical potentials 
H = (/if + /if)/2 and h — (/if — /if)/2. As a first ap- 
proximation, one understands the trapped gas by break- 
ing up the cloud into small pieces, and assuming that at 
each of these pieces is homogeneous and in local equilib- 
rium. Maintaining equilibrium in the presence of energy 
transport requires T is constant, allowing momentum 
transport requires that the pressure obeys VP = — nW, 
where V(r) is the trapping potential, and allowing par- 
ticle transport requires V/i, = — VV^ for j =t,|- This 
hydrodynamic, Thomas-Fermi description allows one to 
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FIG. 1: (Color Online) Schematic phase diagram of a two 
component Fermi gas as a function of (a) Temperature [T] - 
Polarization [P = (n-f — 7i|)/(ri| (b) Temperature [T] - 

chemical potential difference [h = — fj,±)/2], and (c) chem- 
ical potential \p = (/if — /ij.)/2] - chemical potential difference 
[h = (/if — /ijJ/2]. The equation of state sets a relation- 
ship between P, T, h, and /x, so only three of them are needed 
to specify the state. Solid lines: continuous phase transitions; 
dashed lines: discontinuous; these meet at the tricritical point 
[Pt, ht,Tt] or [/it, ht,Tt]. The gray region in (a) maps onto the 
dashed line in (b), and represents a coexistence region. On 
the BCS side of resonance a < 0, and for sufficiently large 
a, Pi=0 and < P2 < 1. When a > decreases in mag- 
nitude, Pi and P2 move to the right, sequentially hitting the 
maximum allowed value P = 1. At unitarity, a = 00, a Wilso- 
nian RG theory[6(| predicts P t = 0.24 and T t /T FA = 0.06. 
Monte-Carlo calculations suggest T C /T F = 0. 152(7) [gj, and 
P 2 = 0.39 HH. 

directly read the structure of the trapped gas from the 
homogeneous phase diagram in figure [He): h is constant 
and \i varies from a large value at the center of the cloud 
to a small value at the edge. The iso-density contours of 
the cloud follow the iso-potcntial contours. 

Observations at Rice are inconsistent with this lo- 
cal density approximation [62| . As already empha- 
sized, the missing element is that the surface energy de- 
scribed in Eq. {T]) distorts the cloud. Microscopically, 
this surface tension arises due to the energy cost of the 
supcrfluid-normal boundary, where the order parameter 
varies rapidly. Depending on the size of the superfluid 
region it is energetically favorable to either shrink this 
boundary to reduce its area or, because the surface en- 
ergy depends on density, shift the boundary to a low 
density region. In a spherical trap this effect changes 
the radius (and density) of the superfluid core. In an 
elongated cloud, one generically expects that the aspect 
ratio of the superfluid region is reduced. Phenomeno- 
logical models based on this principle SH , seem to 



account for the experimental observations. 

By fitting various models to experimental density pro- 
files, Haquc and Stoof [4l[ and DeSilva and Mueller [|o| 
both made an estimate of r\. The quoted values of r\ in 
po| were unintentionally scaled by a factor of fku z / '/xo- 
When corrected for this factor, they found that r\ « 0.6 
for the relatively high temperature data in [l[. Haque 
and Stoof [4l[ found r\ « 4.8 for the lower temperature 
data in (2j. This is consistent with the expectation that 
surface tension should drop as one approaches the tricrit- 
ical point. Here we compare our model with the data in 
0, finding r) w 3. We attribute the slight discrepancy 
with Haquc and Stoof to trap anharmonicities (which wc 
did not include in our calculation). We emphasize that 
?y ~ 3 is more than one order of magnitude larger than 
the value predicted by our microscopic model, r\ « 0.17. 



II. CALCULATION OF SURFACE TENSION 

In this section we present a calculation of the surface 
tension in the BEC-BCS crossover. First we give a very 
crude order of magnitude estimate of the surface tension 
at T = in a unitary gas, then we use a numerical so- 
lution to the BdG equations to obtain a more rigorous 
result at zero temperature. In the deep BEC limit and 
at finite temperature one can use a simpler theory where 
the free energy is expanded in gradients of the order pa- 
rameter. We use this theory to obtain the temperature 
dependece of the surface tension in the unitary limit and 
to test our numerical solution of the BdG equations in 
the BEC limit. 

Related microscopic calculations have been performed 
by Caldas [H| and Imambekov, Bolech, Lukin and Dern- 
ier [H. 



A. Order of magnitude 

Before presenting a detailed microscopic calculation of 
the surface tension we give a simple estimate of its mag- 
nitude at zero temperature in a unitary gas. In the stan- 
dard semi-phenomonological model for surface tension, 
one considers the spatial variation of an order parame- 
ter: in this case the superconducting gap A(r). At the 
first order phase boundary, the free energy has two local 
minima: one with A = Ao and the other with A = 0. 
Surface tension can be attributed to the fact that in the 
boundary between the two bulk phases, A must pass over 
a free energy maximum. 

The energy cost per unit area of the boundary a is 
most simply estimated as the product of the maximum 
height of the free energy barrier per unit volume Sfl and 
the thickness of the domain wall £. The healing length 
£ arises from a competition between the stiffness of the 
order parameter and the height of the energy barrier, and 
at unitarity should be of order the interparticlc spacing. 
Within BCS theory (reviewed in detail below) , the energy 
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barrier is Sfl « 0.2 x n 5 s /3 h 2 /2m. With £ w nj 1/3 
gives 77 w 0.2. As previously quoted, the full solution of 
the Bogoliubov de Gennes equations yield 77 ss 0.17. 



B. Mean Field Theory 

To calculate the properties of the supcrfiuid-nornial 
boundary, we numerically solve the Bogoliubov de 
Gennes equations with appropriate boundary conditions. 
The homogeneous version of these equations is often used 
to describe the zero temperature BEC-BCS crossover in 
ultracold fermions. We emphasize however that near uni- 
tarity, quantitative predictions of this theory should be 
viewed with some skepticism. One could argue that since 
the mean-field theory overestimates the density discon- 
tinuity at the supcrfluid-normal boundary it should also 
overestimate 77. In this section we review the formal the- 
ory, while in the following sections we report results, and 
describe some simplifying approximations. 

Our model consists of a Fermionic system with two dif- 
ferent hyperfme states labeled by a =|, J, in three spatial 
dimensions. The atoms interact via a point interaction. 
The system of N = iVj + iVj_ atoms is then described by 
the Hamiltonian H = J d 3 r(Ho + H- mt ), with kinetic (i?o) 
and interaction (Hi nt ) energy densities 

F (r) = $>J(r)( - |-V 2 - 7v)W(r) 

(4) 

ff int (r) - -^{(r)^(r)^(r)^ T (r) 

where Vv(r) are usual Fermionic field operators. Param- 
eters m, and \x a are the mass and chemical potential of 
the atomic species a respectively Following convention, 
we take | to be the majority species of atoms and use 
variables h = (/i-j- — /ij.)/2 > and fi = (/if + /ij)/2. 
The bare coupling constant U is renormalized by ex- 
pressing it through the s-wave scattering length a s as 
1/U = 1/U R + l/F£ 3 l/2e° with U R = -4irh 2 a s /m, 
e° = h 2 q 2 /2m where V is the system volume [65j . 

Performing a mean-field decoupling of the interaction, 
one writes the Hamiltonian in terms of a gas of Bogoli- 
ubov excitations 16611. 



* T (r) 
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(5) 



Excitations with spin a =|, J, have energies E nM = E n ± 
h, where E n is the positive energy eigenvalue of the BdG 
equations 



-S-A* A(r) 




u„(r) 



(6) 



The mean-field free energy is then [67j, [68| 

|A(r)| 2 



D.MF 
An 



= { £n ~ ^ ~ A «} + / d3r 

n 



log 2 cosh ; 



U 



log 2 cosh 



(7) 



where 6n ~ I 1 



Tr 



V 2 



and (3 = l/k B T. 



The chemical potentials are set by the number equations 
N a = J d 3 rn a (r) with 

n a {r) = (V4(r)^(r)) (8) 
= K(r)\ 2 f(E n , a ) + Mr)| 2 [1 - f{E n _ a )] 

a 

where f(E) = 1/(1 + e~P E ). Self-consistency requires 
that the gap obeys 

A(r) = E/(V> T (r)Vu(r)) 

= ^«n(r)<(r) [1 - f(E nA ) f(E nd )} . ( 9 ) 

n 

It is useful to draw attention to three particular limits 
of these equations. First, if A(r) = one has just a 
non-interacting Fermi gas, with free energy 



n 



(10) 



k.rr 



This equation highlights the most significant flaw of the 
mean-field approach, namely that it yields a noninter- 
acting normal state, overestimating the stability of the 
supcrfluid. 

As a second useful limit, one can consider the case 
where A(r) = Ao is uniform. In that case one can label 
the eigenstates of (J6j) by momentum and one finds the 
standard result 



^hom 



= E^-M-A fc } + F-^ 



(11) 



log 2 cosh 



0(Ek + h) 



log 2 cosh 



P(E k - h) 



where V is the system volume and E^ = 
\J [k 2 /2m — /i) 2 + Aq. In this uniform limit one 
finds that the total particle number N = N-[ + Ni, 
population imbalance AA = Af — Af and gap equation 
become 



5 



0.8r 



0.4 - 



0.1 - 



0.0 - 




1.0 



0.8 



0.6 



0.4 



0.2 



0.0 



▲ Al/fcpa=0.35 



i 

A 



A 
A 

A 
A 
A 




20 



FIG. 2: (Color Online) Order parameter profiles at the interface between normal and superfluid at critical Zeeman field h c . 
Left to right: BCS to BEC side of resonance. Each data point corresponds to a single gridpoint of our real space discretization. 
Insets: normal state T-matrix (pair susceptability) as a function of momentum q at the first order phase transition line h = h c 
corresponding to the same parameters as the BdG calculations. The Fourier transform of T(q) describes the decay of the 
superfluid order parameter into the polarized normal state. The vertical line shows q = k F — kp. 
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(12) 
(13) 

(14) 



A third important case is when A(r) is periodic. Pro- 
totypical examples are the FF state A = Appe %qx or the 
LO state A = Alo cos(qx). In the BCS limit there is 
a small window of stability for such a solution. These 
exotic superfluid states can be described by the BdG ap- 
proach. 

In order to compute the surface tension across the 
BEC-BCS crossover wc numerically solve the BdG equa- 
tions in Eq. ([6]). We find the parameters such that bulk 
normal and superfluid have the same free energy and then 
minimize the free energy functional in Eq. ([7]) with re- 
spect to the order parameter for a domain wall between 
the two phases. A simple way to minimize the free en- 
ergy functional is to iterate the gap equation Eq. (J9j> to 
self consistency. We find that it is more computation- 
ally efficient to directly minimize Eq. ([7]) with respect to 



a discretized representation of A(r). For efficiency we 
calculate the gradient using 



<«W[A(r)] 



SA(r) 

1 7 
U 



2 
U 



[A(r) - A(r) 



(15) 



A(r) = ^ UB (rK(r)[l-/(£„, r )-/(£ n , 



where we used 5E n /SA(r) = 2u n (r)v n (r). 



C. Results (T = 0) 

The domain wall profiles calculated within the BdG ap- 
proach are shown in figure [21 As can be seen, on the BEC 
side of resonance the domain wall is largely featureless, 
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while oscillatory structures develop as one approaches the 
BCS limit. In section Hi D 21 we relate these features to 
the behavior of the T-matrix, and draw the connection 
between these oscillations and the FFLO phase. 

Due to the contribution from the interface, the energy 
found in these calculations exceeds that of the bulk super- 
fluid/normal gas. From this excess energy we extract the 
dimcnsionlcss parameter ?/; our results are summarized 
in Fig.[3l At unitarity we find 77 ~ 0.17. The surface ten- 
sion drops as one approaches the BCS side of resonance. 
It grows to a maximum of 77 ~ 0.25 at 1/fc/a w 0.4, then 
falls as one proceeds towards the BEC limit. 
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FIG. 3: (Color Online) Dimensionless surface tension r\ = 
2h~ 2 mn^ 4 ^ a a as a function of (fc^a) 1 at T = 0. When 
(fc^a)" 1 > 1.01 the superfluid state is partially polarized. 
Triangles: calculation using the full BdG equations as de- 
scribed in III CI circles: gradient expansion approximation to 
this solution from III Dl The lines are a guide to the eye. 
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D. Gradient expansion 



Even in the variational form discussed in the previous 
section, solving the BdG equations is numerically chal- 
lenging - especially at finite temperature. To explore how 
the surface tension varies with temperature we consider 
a Landau- Ginzburg expansion, 



/ d 3 r [ai |A(r)f + a 2 |VA(r)|* + a 3 |V 2 A(r)|" + 

7l |A(r)| 4 + 72 |A(r)| 6 + --- 
gradient terms of order higher than A 2 

(16) 



The sum of the terms a±, and 71, 72,'' '7n yields the 
homogeneous free energy in Eq. pip , and we include all 
of them by numerically calculating f^hom- Additionally, 
in this section we calculate the gradient terms which are 
quadratic in A, namely a 2 , 0:3, ■ ■ ■ . We will neglect gra- 
dient terms higher than 0(A 2 ), an approximation well- 
justified near the tricritical point. By comparing to our 
solutions of the BdG equations, we find that this ap- 
proximation introduces about a factor of 2 error at zero 
temperature. For the terms which are quadratic in A, we 
go to all orders in the gradient. At low T the first term, 
«2j becomes negative and higher order terms stabilize the 
system. 

This expansion is readily developed from the path inte- 
gral representation of the partition function, after decou- 
pling the interaction term with a Hubbard-Stratonovich 
transformation 



Su 



rP rP 
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(17) 
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where the expectation (e Sl ) is evaluated in the free partition function in terms of the pair susceptibility x(<?) : 
fermion ensemble. Evaluating this expectation yields the 

v(«)|A„| 2 +0[A 4 ] 



*(?) = £ 



1 - /(ej +fcT ) - /(ej- fci ) _ J_ 

2e k 



Cs.-i.kt + ea 



fcT I" fc f-fcl 



(18) 
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The function x(o) 1S °f general importance in the 
many-body problem of the BEC-BCS crossover of spin 
polarized fermions, such as T-Matrix approximation 
schemes [65j and the finite-momentum pairing instab ility 
of the polarized normal phase to the FFLO phase |22l.l69j. 
It is the static limit of the more general two-particle prop- 
agator x(q, w n ) 



x(q,M/„) = ^2 



l-/(e§+k, T J -f[e%-k,l 




At zero temperature the integrals evaluate to 

f\ 1 ft I fl I fl 

2^x(q,^) = 2mC-2(l + r)- J+ ^ J -[ J+ J - (20) 



<7 - 2 + 2c< 
9 + 2 + 2ctC 
<2 - 2r + 2aC 
q + 2r + 2a( 



fl = 




/i - 





crgC log 



where x = {^/mVk F )x, q = H/k F , r = k l F /k F , 



mv/(k F ) 2 ,C, = v + (1 + r 2 )/2 - q 2 /4 and k F ' 



1 • I 



i/2m(/J ± h), and for brevity we have set h = 1. When 

A*t 1 — 0, one takes fct = 0. 

The normal state becomes locally unstable towards fi- 
nite q pairing when the coefficient of |A g | 2 in Eq. (fTS)) 
becomes negative: 



x(q,0)--=0. 



(21) 



Above the tricritical point this Thouless criterion locates 
the second order phase boundary, while below the tricrit- 
ical point it yields the spinodal of the first order phase 
transition. 



1. Mean Field approximation 

At this point we neglect fluctuations in A, finding that 
to quadratic order 



Z = e-^ a cxp ^-/^T-^IAJ 2 ^ 



(22) 



where the T-matrix obeys T^ 1 (q) = l/U — x(l)- Thus 
explicit expressions for a±, a^i a 3>' ' ■ are found by ex- 
panding x(q) i n powers of q. 

Using this expression for the quadratic terms in (|16l) . 
we get 



fi[A] 



,3 ^hom(A(r)) 

ir - 



-J2(x(0)- X (q))\A q \ 2 , (23) 



where A(r) = V 1 S g e" 1 ' r A g . The q = term is sub- 
tracted to eliminate double-counting of the ai term. 



We calculate the surface energy from (|23|) by employing 
the ansatz 



A(z) = A (erf(4z/IF dw ) + l)/2 



(24) 



and minimizing Q with respect to Wd w , which is a mea- 
sure of the width of the domain wall. Technical details 
are describe in appendix [B] Figure [3] shows the temper- 
ature dependance of the domain wall width at unitarity. 
As can be seen, the domain wall size diverges at the tri- 
critical point. 




FIG. 4: Width (Wdw) of the domain wall as function of inverse 
temperature (/3) at unitarity, with parameters measured in 
terms of the chemical potential /i and k% = 37r 2 n s . Top scale 
is nonlinear. The domain wall width diverges at the tricritical 
point around f3(i ~ 2.0. 

The surface tension ?y, extracted from the free energy 
evaluated at the optimal W^ m , is shown in figure O We 
also show what this surface tension would be if we expand 
x(q) to either quadratic or quartic order in q. This cor- 
responds to truncating (fT6|) at ct2 or While these 
latter approximations work well around the tricritical 
point, they do not correctly describe the low tempera- 
ture physics: both «2 and change sign at low temper- 
ature, and without the influence of higher order terms, 
the normal state becomes unstable to a FFLO state, and 
the surface tension vanishes. Using the full x, we find 
that as T — > 0, the dimcnsionlcss surface tension be- 
comes rj ~ 0.1. The discrepency between this result and 
the full BdG equations can be attributed to the neglect 
of gradient terms which are higher order in A. 

Comparing these gradient expansion results to the full 
solution to the BdG equations at T — 0, we find that 
the agreement becomes better in the deep BEC limit 
(see Fig. [3]). In this limit, x is weu approximated by 
a parabola, and fihom can be truncated at quartic or- 
der in A. Minimizing the resulting free energy results in 
a Gross-Pitaveskii equation [7(j. Recently, Sheehy 71 1 
has discussed the role of quantum fluctuations near the 
tricritical point. 
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FIG. 5: Dimensionless surface tension n = 2h~"mn a ~'~o as 
a function of inverse temperature fi, keeping different num- 
bers of terms in the gradient expansion. Even with the full 
X, we include terms up to C(A ) in the free energy. As tem- 
peratures goes to zero, the solid line suggests n ~ 0.1. This 
should be compared with the full solution of the Bogoliubov- 
de-Gennes equations which give rj — 0.17. The kink near 
fill ~ 6 is an artifact of neglecting gradient terms which are 
beyond quadratic order in A. 

2. Proximity effect 

On the normal side of the domain wall, where A(z) 
is small, the free energy can be expanded to quadratic 
order in A. To find the asymptotic behavior of A in the 
domain wall one minimizes the quadratic approximation, 
fl « Eq T_1 (q)l A q| 2 witn tlic constraint that A(z = 
0) = A is nonzero. By symmetry, A q 
finds 



A qz , and one 



T- 1 (q z )A qz -\ = 0, 



(25) 



where A is a Lagrange multiplier. Consequently, for large 
z, the order parameter is proportional to the one dimen- 
sional Fourier transform of the T-matrix 



A(z) cx J ^e*T{q). 



(26) 



This result is confirmed in figure [2l where insets show 
the behavior of T(q). For example, on the BEC side of 
resonance the T-matrix is peaked at q = 0, yielding a 
monotonically decaying order parameter. On the BCS 



side, the T-matrix is peaked at finite q ss k F 



kp, and 



one observes oscillations of A(z) with this wave-vector. 

One can make a simple analogy, noting that the q- 
dependance of the T-matrix is reminiscent of the fre- 
quency dependance of a driven damped harmonic oscil- 
lator. The spatial dependance of the order parameter 
domain wall would then be analogous to the temporal 
dependance of the oscillator's position once the drive is 
turned off. The BEC/BCS sides of resonance then map 
on to overdamped/underdamped oscillators. 



III. EFFECT OF SURFACE TENSION ON 
DENSITY PROFILES 

Having explored the microscopic theory of the domain 
wall separating the superfluid and normal regions, we 
now investigate how the surface energy in this bound- 
ary affects the shape of a trapped gas. We will assume 
that the zero temperature population imbalanced atomic 
system is phase separated into two regions: a central su- 
perfluid core surrounded by a normal shell. We will take 
the normal state to be fully polarized (with = 0) and 
the superfluid state to be fully paired (with = njj. As 
discussed in the introduction, this is an approximation. 
Remarkably, the experiments at Rice university [l], are 
largely consistent with this ansatz, which was first in- 
troduced by Chevy [z3] . As seen in figures [5][71 in these 
elongated clouds there is no experimental evidence of a 
partially polarized normal region between the fully paired 
superfluid and fully polarized normal regions. The ab- 
sence of this phase is significant: a partially polarized 
normal state is seen in experiments at MIT 0, 0, H[ and 
in more sophisticated theoretical calculations of the bulk 
phase diagram [2^ . IBH . [731 ]. The unexplained behavior 
seen at Rice is presumably related to the small numbers 
of particles and the high aspect ratio of the cloud, but 
other considerations, such as the kinetics of spin trans- 
port, may also be important (3rl l38l. [7l. |75| . 

We will restrict our discussion to unitarity, where 
physics is universal and the superfluid and surface en- 
ergy densities between the superfluid and normal regions 
have simple forms. The equation of state of the central 
superfluid shell is given at T = by eq. ©, while the 
outer fully polarized normal shell obeys 



?i„ fr) = rit = — ^ 
6tt 2 



2m/i„(r) 



K 2 



3/2 



(27) 



The free energy densities of the bulk phases f s 
— J n Stn diJ, can be written as 



fsA r ) = - 



15tH 



2m 



3/2 



(28) 



where ( s = 1/(1 + /3) 3 / 2 , Q n = 1/2. Then we calcu- 
late the total bulk energies CIs,n = J s / n d 3 'rf s , n [n(r),h} 
by integrating the bulk energy densities over the super- 
fluid/normal regions. As previousy introduced, we model 
the surface energy density a[fi(r),h] — i](h 2 /2m)ni^ 3 , 
where rj is the dimensionless constant calculated in the 
previous section. We calculate the total surface energy 
Edw = J d 2 ra[/j,(r), h] by integrating the surface en- 
ergy density over the superfluid-normal boundary. Away 
from the superfluid-normal boundary, we assume that 
the system is locally homogenous and the external har- 
monic trapping potential Vt, ap (r) = b±p 2 + b z z 2 = 
mu> 2 (X 2 p 2 + z 2 )/2 is treated in the LDA by introducing 
a local chemical potential /i(r) = fio — Vt rap (r). Given 
that the experimental traps arc formed by focussed laser 
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FIG. 6: Experimental two-dimensional column densities (black denotes high density) for P = 0.38 with theoretically calculated 
boundaries for different surface tensions rj (fixing the number of particles to be constant). Top: majority atoms 7Vj; Bottom: 
minority atoms Afj.. The dotted line is the ellipse with semi-major and semi-minor axes Ztf and -Rtf respectively, while 
the solid line is the superfluid-normal boundary in the presence of surface tension. As rj is increased, the superfluid-normal 
boundary deforms from an elliptical isopotential surface, but the boundary becomes increasingly insensitive to surface tension 
with increasing rj. N c = 15 Fourier components were chosen for equation (|30f) . Data corresponds to Fig. 1(c) in Ref. 0], 
used with permission. Data outside of an elliptical aperture has been excluded. This truncation of the data leads to a slight 
discrepancy in P compared to the value quoted in [jj]. Each panel is 1.4mmx0.06mm, and shows the true aspect ratio of the 
cloud. 



beams, describing it via a harmonic potential should be 
viewed as an approximation. 



A. Calculation of boundary 

We make a completely general ansatz for the domain 
wall, only assuming rotational symmetry about the long 
axis of the trap. We parameterize the boundary in terms 
of coordinates / and 9, which are related to the cylindri- 
cal coordinates p and z by 

p(ej) = r tf j cos e 

z(6J) = Z TF f sm9 1 ' 



where Rtf = \/po/b±, Ztf = \J Ho/b z - The bound- 
ary is described by the function / = F(9). As shown in 
appendix[Cl the two-dimensional integrals for the free en- 
ergy can then be simplified to one dimensional integrals, 
which can be performed numerically. 

The optimal shape is found by minimizing the free en- 
ergy functional £It = &s + + -E-xw on the space of 
functions F(9) at fixed JVj and N^. The constraints are 
imposed using Lagrange multipliers. 



We expand F{9) as 

00 

F{9) = an cos(2?7.(9) (30) 

71=0 

which is compatible with the boundary conditions im- 
posed by the symmetry of the problem, F'(0) = 
F'{n/2) = 0. We truncated this series at a finite number 
of Fourier components N c and numerically minimized f2y 
with respect to ao, 01, . . . , a/v c - We find that we need to 
include more terms in this series when 77 is larger, but for 
all values of 77, the profiles become insensitive to N c for 
N c > 15. 

In figure [6] we plot the boundary F{&) that minimizes 
VLt for different values of r/. The boundary becomes 
almost insensitive to 77 for high surface tension. This 
behavior has two sources: (i) For large 77 the ends be- 
come increasingly flat, so surface tension plays an increas- 
ingly insignificant role, (ii) the edges along the minor 
axis touch the edge of the majority cloud, at which point 
the superfluid-normal boundary changes to a superfluid- 
vacuum boundary and surface tension ceases to be im- 
portant. Due to this "saturation" of the boundary shape 
with high 7], and the difficulty of defining the bound- 
ary from noisy 2-D data, we find it convenient to fol- 



10 




FIG. 7: Axial densities. Symbols: experimental one-dimensional Li spin densities and density differences for P — 0.39 
(JV t = 155,000, N t = 68,500) (left column) and P = 0.63 (iV t = 123, 600, JV| = 28,000) (right column), from Ref. 0, with 
permission. Lines: theoretical curves for r\ = 2.83, taking a cigar shaped harmonic trap with small oscillation frequencies 
lj z = (2-7r)7.2Hz and ui r = (27r)325 Hz. Oscillations in the density difference within the superfluid region are artifacts of our 
ansatz (|30|) . To minimize noise, only experimental data inside an elliptical window was considered (see text). This aperture is 
visible in figure [6] 



low references 5l[ and find r] by fitting our the- 
oretical model to the 1-D axial densities, defined by 

n \ a \( z ) = § dx dy n(x,y, z). As illustrated in Fig. [6j we 
improve signal to noise by excluding data outside of an 
elliptical window (77[. We find that r\ ~ 3 gives an ax- 
ial density difference profile most closely matching the 
experimental density from Ref. Q for P = 0.38 and 
P = 0.63. As seen in figure [71 the overall quality of the 
fit is quite good. There are however distinct differences 
between the predictions of the model and the observed 
profiles. These can largely be attributed to trap anhar- 
monicities whose modelling is not reported here. 

We also believe that the Sp dependance of rj may be 
important for capturing the exact shape of the domain 
wall. Generically one would expect that this dependance 
would reduce 77 at the ends of the boundary, increasing 
the curvature of the end-caps and making a smoother 
axial density. This effect would also lead to an apparent 
polarization and number dependance of ?/. Finally, we 



found some sensitivity to how we treat the background 
in each image. For example, if we fit the axial density 
difference at P = 0.6 without windowing the data, we 
find that rj = 1 provides a better fit. This sensitivity 
can be attributed to structure in the background which 
persist throughout the image, even far from the cloud. 

Since they are based upon identical models (just using 
different ansatz's for the boundary shape), the quality 
of our fits are very similar to the ones found by Haque 
and Stoof when investigating a large number of similar 
profiles fill ]. Converting to our units, Haque and Stoof 
found 77 = 4.8 ± 1.2. Their result is slightly higher than 
ours. We attribute this difference to differences in fitting 
procedures (such as windowing the data) and to modeling 
of the trap. Haque and Stoof used a more sophisticated 
Gaussian model for the trap, while wc assumed it was 
harmonic. 

The authors of Ref. 0] from the MIT experiment find 
no visible distortion of the superfluid cloud and quote 
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an uncertainty of about 2% for this null observation. We 
studied how surface tension would affect this experiment. 
As evident from figure [H a distortion of less than 2% im- 
plies rj < 1. Thus their null-observation is consistent with 
a value of surface tension on the same order as what we 
find in our microscopic calculation. Figure [8] shows that 
the surface tension needed to fit 0, rj ~ 3, would cause a 
distortion more than 10% in the MIT experiment (as was 
already pointed out in [f|), well above their threshold of 
detection. 



a) 

a 











































































































For MIT experiment: 

iV t =6-10 fi , P=Q.U, A=(U5 











i.O 0.5 1.0 1.5 2.0 2.5 3.0 3.5 

V 



FIG. 8: (Color Online) Distortion of superfluid core aspect 
ratio (= 1 — F(n /2) / F(0)) in % as a function of the dimen- 
sionless surface tension rj for parameters of Ref [f|, where A 
is the aspect ratio of the harmonic trap. 



IV. SUMMARY AND CONCLUSIONS 

We presented a microscopic calculation of surface ten- 
sion at a superfluid-polarized boundary in a polarized 
two-component Fermi gas, finding a value of the dimcn- 
sionless surface tension r\ which is consistent with an e- 
expansion [76j . We analyze both the interaction depen- 
dance and temperature dependance of r\. We argue that r/ 
vanishes in both the BEC and BCS limits, with a peak on 
the BEC side of resonance near unitarity. Furthermore, 
rj decreases with temperature, vanishing at the tricritical 
point. 

We compared our microscopic calculation of rj to ex- 
perimental estimates extracted by fitting a phenomeno- 
logical model to axial density data. We find that our 
microscopic predictions of r/ = 0.17 at unitarity are an 
order of magnitude smaller than our best fits to the ex- 
perimental data from Rice (77 ~ 3). On the other 
hand, phenomenological modeling of the the MIT exper- 
iment [5j bound 77 to be somewhere between zero and a 
few times larger than our microscopic predictions. We 
think that additional theoretical insight, as well as more 
experimental data, would be needed to resolve these dis- 
crepancy. 



Despite the differences in the magnitude of 77, we find 
that the experimental observations at Rice Q show all of 
the appropriate hallmarks of surface tension. Using our 
phenomenological model, we presented a short study of 
how surface tension distorts the shape of the superfluid 
core in a harmonically trapped cloud of atoms. We repro- 
duce both the observed double peak structure [l[ in the 
axial density difference profile and the distorted shape of 
the supcrfluid/normal boundary Q. 

We find that with increasing surface tension (fig- 
ure [6]), the superfluid/normal boundary attains a "lim- 
iting" shape significantly different from the isopotcntial 
contours that are predicted by the Thomas Fermi ap- 
proximation. As the temperature increases, the system 
moves closer to the tricritical point in figure [T] and as a 
result the surface tension decreases and disappears at the 
tricritical point in figure [5] As a result, surface tension- 
induced distortion of the superfluid core must be ab- 
sent if the temperature of the atomic trap is maintained 
above the tricritical point; indeed, this behavior is ob- 
served in We would argue that thermal effects are 
not responsible for the differences between experiments 
at Rice [2, @] and MIT [1,0, H- The strongest evidence 
that temperature is not the issue, is the excellent agree- 
ment between zero temperature Monte-Carlo results [59| 
and the MIT experiments Q . 

Sensarma et al (36[, and more recently Tezuka and 
Uedajll], attempt to understand the deformation of the 
superfluid core by studying a microscopic model of the 
entire atomic cloud. They solved the BdG equations for 
a relatively small number of particles (a few thousand) in 
an axially symmetric system. While these calculations do 
give more insight into the properties of these systems, it 
is difficult to extract quantitative information from them. 
In particular, the small particle numbers lead to a much 
larger surface area to volume ratio than in experiments, 
and artificially amplifies the role of surface tension. Wc 
believe that those calculations do not "explain" the in- 
consistencies which we observe between the magnitude 
of surface tension in our microscopic BdG calculations 
and our phenomenological modeling of the experimental 
data. It is possible that larger simulations of that sort 
might be more useful in this regard. 
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APPENDIX A: CUTOFF DEPENDANCE OF THE 
BDG CALCULATIONS 

In order to compute the surface tension across the 
BEC-BCS crossover we used a numerical calculation 
based on the BdG equations (|6|). Here we give details 
about our numerical approach, showing that we have 
used sufficiently large cutoffs to produce unbiased results. 

We find parameters such that bulk normal and super- 
fluid have the same free energy and then minimize the 
functional ([7]) with respect to the order parameter for a 
configuration with a domain wall between the two phases. 
The excess free energy of this configuration is attributed 
to the domain walls, and allows us to extract a surface 
tension. 

We assume that the order parameter varies only along 
the z-direction. In the simulation we discretize this one 
dimensional space on a uniform grid with TV grid points 
and approximate the gradients through a fourth order 
finite difference matrix. We find it convenient to use pe- 
riodic boundary conditions in the z-dircction, simulating 
two domain walls. We have verified that the interaction 
between the two walls is negligable. 

Translational invariance in the directions perpendicu- 
lar to the interface(i.e. x- and y-directions) implies that 
the u and v'a are of the form 



u(r) 



ikj_rj_ 



u n (z) 



V(Y) 



-v n (z) (Al) 



One has to solve a 27V x 27V matrix eigenvalue problem 
for each k± = |kj_|. The coupling constant U depends on 
the UV-cutoff and was fixed by rcnormalizing Qmf{^-o) 
through the result from a direct calculation (fTTj) where 
a uniform order parameter was assumed. We found that 
when seeded with an order parameter profile with two 
domain walls the minimization algorithm converges to 
a local minimum with two domain walls. This solution 
correctly obeys the self-consistent gap equation. 

We systematically increased TV to check the conver- 
gence of the order parameter profiles and surface tension 
(see Fig. [5]). In the deep BEC limit, the gradient expan- 
sion becomes a good approximation to the BdG equa- 
tions, enabling us to check our BdG calculation. In that 
limit we found excellent agreement between these theo- 
ries, giving us confidence in the accuracy of our numerical 
methods. The largest number of grid points TV and trans- 
verse modes TV& we used were TV = 400 and TV^ = 400; 
these were, for example, used to obtain Fig[3J We feel 
that any residual errors from our finite gridpoints/cutoffs 
are much smaller than the errors introduced through the 
mean field approximation. 



APPENDIX B: DETAILS OF GRADIENT 
EXPANSION 

This appendix contains technical details of the gradi- 
ent expansion used for finding the domain wall energy at 
finite temperature. 
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FIG. 9: (Color Online) Top: Representative order parameter 
profiles for different cutoffs computed using the BdG equa- 
tions at I/Zcfh = 0.05. For better visibility a line connecting 
the data points is displayed. Bottom: Dimensionless surface 
tension constant r\ for different cutoffs as a function of 1/kFa. 
(Color online) 



We calculate the uniform gap (Aq) by solving 



dy 



1 - 



sinh E„ 



y 



cosh_E.y + cosh (3h E y 



= 



(Bl) 



where E y = ■J (y 2 — /3/i) 2 + (/3A) 2 . We numerically per- 
form the integral up to a cutoff y c , then analytically ap- 
proximate the remainder by expanding the integrand to 
sixth order in 1/y. We calculate the free energy of the 
solution from (J7J), again breaking the integral into two 
pieces, and analytically integrating the tail. We then 
find the parameters for which the normal and supcrfluid 
states have equal energy. 

We write {2k B T/V)x(q) = Ax{q^/fi/2m), with A = 
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l/(27r 2 )(2m//3) 3 / 2 to express 

kC(k, q) - 2k 2 - \q 2 



x(q) 



dk- 



k 2 ~[3p 

C(k,q) = -\oz{Z\Z i _jz\z\ r )-2k 

Z° — l _|_ e (fe-sg/2) 2 -/3/i<r 



(B2) 

(B3) 
(B4) 



After tabulating this integral, calculating the energy of a 
domain wall is straightforward within the error function 
ansatz of Eq. ([24]) . Again, we split all integrals into two 
parts, integrating the tails analytically. 

To extract r\ we also calculate the bulk supcrfluid den- 
sity, 



1 (2rn\ 3/2 



dy y 2 x 



sinh E„ 



E y J cosh E y + cosh f3h 



(B5) 



using A = A . 



APPENDIX C: EVALUATION OF 
PHENOMENOLOGICAL FREE ENERGY 

This appendix gives the analytic expressions used to 
calculate the free energy of an arbitrary domain wall. 
We parameterize the boundary by / = F(9), in terms 
of which the coordinates of the boundary are p(9, f) = 
RtfJ cosf and z(6, f) = ZtfJ sin#. The surface energy 
is 



Edw — A c iw 
= 2A d 



1 - 



tt/2 



p 



E>2 y2 

JXrpp TF 



d9F(9) COS 9 

x [F'{6) cos 9 - F{9) sin 9} X [1 - F{9) 2 ] 2 

\ ( Z TF \ 2 ( F>(6)sm9 + F{9) cos U 
X V + \Rtf) \F'{0) cos9-F(9) sin 9 



(CI) 



where we define the coefficient 
Adw = It^JzRtf^tf 



"2m" 


3/2 


h 2 





3/2 



(l+P) 2 ^ 2 )^ 3 



(C2) 



We write the free energy of the supcrfluid core as 

„2 ,2 \ 5/2 



fig = A s pdpdz 1 



p2 72 
^TF TF 



,tt/2 pF(6) 

2A S d9cos9 d// 2 (W 2 ) 5/2 
Jo Jo 

fiv/2 



2A S / d9Gx[F (9)} cost 
Jo 



A. 



( s R tf Ztf 



2m 



-1 3/2 



ft 2 



A S/ 2 

w 

157T 



(C3) 
(C4) 



d(x) = — x\/l - x 2 (-15 + 118x 2 - 136.T 4 + 48a 
384 



+15 sin 1 {x) 



(C5) 



Similarly, the free energy of the fully polarized normal 
shell, Q N = A n J n pdpdz{\ + 7 - p 2 /R 2 TF - z 2 /Z 2 F ) 5 / 2 , 
is: 



tt N = 2A„(l + 7 ) 4 

A n = —C, n R TF ZTF 



5tt 


.tt/2 


r f{9) - 






/ dOGi 


cos 9 


256 ~ 


Jo 


Lvi+7. 



2m 



1 3/2 



h 2 



A 5/2 

15tt 



(C6) 



The total number of atoms in the two phases are given 
by 



/•ir/2 

N S = 2B S d9 cos 9G 2 [F \9)] 
Jo 

N n = 2B n (l + 7 ) 3 



7r/2 

tt/32- / d9 cos 9G 2 




F(9) 



v/r+7 

(C7) 



where 



B a n — Cs.n r, 

Sir 



2mp Q 



r, 2 



3/2 



R\ f Ztf 



G 2 (x) 



1 

48 



Vl -x 2 (-3 + 14x 2 - 8x 4 ) + 3sin _1 (x) 

(C8) 



Thus both the free energy, and the constraints of fixed 
N and P reduce to one dimensional integrals. 
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